Earth to Mars trajectory and departure conditions
The following examples show the use of CelestLab functions for basic interplanetary trajectory calculations. The trajectories are computed using Keplerian propagation (central force only) and patched conics approximation.
The objective is to determine an Earth to Mars trajectory and all departure conditions from a low Earth circular orbit.
The Hypotheses are:
- Departure date from Earth: 1 March 2016
- Arrival date at Mars: 1 August 2016
- Altitude of Earth injection orbit: 500 km
Note: First execute the code below to initialize some variables and utility functions.
// ----------------------------------------------- // Data // ----------------------------------------------- // Gravitational parameter for Sun MUSUN = CL_dataGet("body.Sun.mu"); // Gravitational parameter for Earth = default value MU = CL_dataGet("body.Earth.mu"); // Dates of departure from Earth / arrival at Mars (TREF) cjd_dep = CL_dat_cal2cjd(2016, 3, 1); cjd_arr = CL_dat_cal2cjd(2016, 8, 1); // Semi-major axis of (circular) injection orbit sma_ini = 500.e3 + CL_dataGet("eqRad"); // ----------------------------------------------- // Utility functions // ----------------------------------------------- // Plot trajectory in ecliptic (2D) function Plot_traj(cjd_dep, cjd_arr, pos_dep, vel_dep) kep_dep = CL_oe_car2kep(pos_dep, vel_dep, mu = MUSUN); cjd = linspace(cjd_dep, cjd_arr, 100); kep = CL_ex_kepler(cjd_dep, kep_dep, cjd, mu = MUSUN); pos = CL_oe_kep2car(kep); AU = CL_dataGet("au"); scf(); plot(pos(1,:)/AU, pos(2,:)/AU, "thickness", 2); xtitle("Trajectory in ecliptic (EOD)", "X (AU)", "Y (AU)"); CL_g_stdaxes(); a = gca(); a.isoview = "on"; endfunction // Plot some orbital elements (returns figure id + bounds for second plot function) function [fig, bounds]=Plot_orbElem(cir) inc = cir(4,:); gom = CL_unMod(cir(5,:), 2*%pi); pso = CL_unMod(cir(6,:), 2*%pi); bounds = [zeros(3,2); 0, %pi; min(gom), max(gom); min(pso), max(pso)]; fig = scf(); subplot(2,1,1); plot(CL_rad2deg(gom), CL_rad2deg(inc), "thickness", 2); xtitle("X = RAAN (deg) - Y = Inclination (deg)"); CL_g_stdaxes(); subplot(2,1,2); plot(CL_rad2deg(gom), CL_rad2deg(pso), "thickness", 2); xtitle("X = RAAN (deg) - Y = Argument of latitude (deg)"); CL_g_stdaxes(); endfunction // Plot additional results function Plot_newOrbElem(fig, bounds, cir) inc = cir(4,:); gom = CL_rMod(cir(5,:), mean(bounds(5,:))-%pi, mean(bounds(5,:))+%pi); pso = CL_rMod(cir(6,:), mean(bounds(6,:))-%pi, mean(bounds(6,:))+%pi); scf(fig); subplot(2,1,1); plot(CL_rad2deg(gom), CL_rad2deg(inc), "ro", "thickness", 2); subplot(2,1,2); plot(CL_rad2deg(gom), CL_rad2deg(pso), "ro", "thickness", 2); endfunction |  |  |
1) Interplanetary transfer
Determine Earth and Mars positions at departure and arrival dates and compute the velocity vectors using Lambert method. All calculations are done in EOD frame supposed to be inertial.
// ----------------------------------------------- // Transfer from Earth to Mars (EOD) // ----------------------------------------------- // Positions of planets from Sun ~ Solar system barycenter // Frame = EOD, center = Sun [posbody_dep, velbody_dep] = CL_eph_planet("Earth", cjd_dep); [posbody_arr, velbody_arr] = CL_eph_planet("Mars", cjd_arr); // Transfer: computes heliocentric velocity vectors at arrival and departure duration = (cjd_arr - cjd_dep) * 86400; // seconds [vel_dep, vel_arr] = CL_man_lambert(posbody_dep, posbody_arr, duration, mu=MUSUN); // Departure hyperbolic excess velocity (EOD) vinf_dep = vel_dep - velbody_dep; // Plot trajectory Plot_traj(cjd_dep, cjd_arr, posbody_dep, vel_dep); | | |
2) Departure conditions
Determine all possible injection orbits (inclination, RAAN, argument of latitude) where the (tangential) injection maneuver has to be performed. The orbital elements are obtained at the perigee of the hyperbolic orbit.
// ----------------------------------------------- // Departure conditions on hyperbolic orbit (ECI) // ----------------------------------------------- // Excess velocity in ECI frame vinf = CL_fr_convert("EOD", "ECI", cjd_dep, vinf_dep); // Generate all possible impact vectors (origin = Earth center) // pinf = circle perpendicular to vinf of radius d = "impact" distance d = CL_ip_flybyParams("vinf", CL_norm(vinf), "rp", sma_ini, "dinf"); M = CL_rot_defFrameVec(vinf, vinf, 3, 3); alpha = linspace(0, 2*%pi, 201); pinf = d * M' * [cos(alpha); sin(alpha); zeros(alpha)]; // All position and velocity vectors at periapsis (hyperbola) [pos_hyp, vel_hyp] = CL_ip_flybyVectors("pvinfd", pinf, vinf, ["posp", "velp"]); // ----------------------------------------------- // Departure conditions on circular orbit (ECI) // ----------------------------------------------- // Velocity on circular orbit + orbital elements v_ini = sqrt(MU / sma_ini); pos_cir = pos_hyp; vel_cir = CL_unitVector(vel_hyp) * v_ini; DV = vel_hyp - vel_cir; // DV, not used cir = CL_oe_car2cir(pos_cir, vel_cir); // Plot results [fig, bounds] = Plot_orbElem(cir); | | |
3) Direct method (knowing the inclination of the injection orbit)
First determine the orbit plane (i.e. RAAN) then the argument of latitude.
// ----------------------------------------------- // Direct method // ----------------------------------------------- // Imposed inclination inc = CL_deg2rad(50); // -- Step 1: look for RAAN such that orbit plane contains vinj // W0: angular momentum when RAAN = 0 W0 = [0; -sin(inc); cos(inc)]; [q1, gom1, q2, gom2, Inok] = CL_rot_defRot1Ax(W0, [0;0;1], vinf, %pi/2, numsol=2); // Note: if Inok <> [] => no solution (or only one) if (Inok <> []) q1 = CL_rot_defQuat(ones(4,1) * %nan); q2 = CL_rot_defQuat(ones(4,1) * %nan); end // -- Step 2: look for PSO (= argument of latitude) // PSO computed using Euler angles: axes = [3,1,3], angles = gom, inc, phi // The 3rd rotation angle phi = pso + %pi/2 + turnang/2: // pso => angle from ascending node to injection position // %pi/2 => angle between radius and velocity vectors (circular orbit) // turnang/2 => cruise on hyperbola from perigee to infinity (velocity // vinf) W = CL_rot_rotVect([q1, q2], W0); q = CL_rot_defRotVec([1;0;0], [0;0;1], vinf, W); angles = CL_rot_quat2angles(q, [3,1,3]); turnang = CL_ip_flybyParams("vinf", CL_norm(vinf), "rp", sma_ini, "turnang"); // Solutions (circular elements) newcir = [[sma_ini; 0; 0] * ones(1,2); // sma, ex, ey angles(2,:); // Inclination = [inc, inc] angles(1,:); // RAAN = [gom1, gom2] angles(3,:) - %pi/2 - turnang/2]; // PSO // Plot additional results in already created figure Plot_newOrbElem(fig, bounds, newcir); | | |